home *** CD-ROM | disk | FTP | other *** search
/ The Very Best of Atari Inside / The Very Best of Atari Inside 1.iso / mint / mntlb20 / lib / atof.c < prev    next >
C/C++ Source or Header  |  1992-05-17  |  20KB  |  848 lines
  1. /*
  2.  *     double strtod (str, endptr);
  3.  *     const char *str;
  4.  *     char **endptr;
  5.  *        if !NULL, on return, points to char in str where conv. stopped
  6.  *
  7.  *     double atof (str)
  8.  *     const char *str;
  9.  *
  10.  * recognizes:
  11.  [spaces] [sign] digits [ [.] [ [digits] [ [e|E|d|D] [space|sign] [int][F|L]]]]
  12.  *
  13.  * returns:
  14.  *    the number
  15.  *        on overflow: HUGE_VAL and errno = ERANGE
  16.  *        on underflow: -HUGE_VAL and errno = ERANGE
  17.  *
  18.  * bugs:
  19.  *   naive over/underflow detection
  20.  *
  21.  *    ++jrb bammi@dsrgsun.ces.cwru.edu
  22.  *
  23.  * 07/29/89:
  24.  *    ok, before you beat me over the head, here is a hopefully
  25.  *    much improved one, with proper regard for rounding/precision etc
  26.  *    (had to read up on everything i did'nt want to know in K. V2)
  27.  *    the old naive coding is at the end bracketed by ifdef __OLD__.
  28.  *    modeled after peter housels posting on comp.os.minix.
  29.  *    thanks peter!
  30.  *
  31.  *    mjr: 68881 version added
  32.  */
  33.  
  34. #if !defined (__M68881__) && !defined (sfp004)
  35.  
  36. #if !(defined(unix) || defined(minix))
  37. #include <stddef.h>
  38. #include <stdlib.h>
  39. #include <float.h>
  40. #endif
  41. #include <errno.h>
  42. #include <assert.h>
  43. #include <math.h>
  44.  
  45. #ifdef minix
  46. #include "lib.h"
  47. #endif
  48.  
  49. #ifndef _COMPILER_H
  50. #include <compiler.h>
  51. #endif
  52.  
  53. extern int errno;
  54. #if (defined(unix) || defined(minix))
  55. #define NULL     ((void *)0)
  56. #endif
  57.  
  58. #define Ise(c)        ((c == 'e') || (c == 'E') || (c == 'd') || (c == 'D'))
  59. #define Isdigit(c)    ((c <= '9') && (c >= '0'))
  60. #define Isspace(c)    ((c == ' ') || (c == '\t'))
  61. #define Issign(c)    ((c == '-') || (c == '+'))
  62. #define IsValidTrail(c) ((c == 'F') || (c == 'L'))
  63. #define Val(c)        ((c - '0'))
  64.  
  65. #define MAXDOUBLE    DBL_MAX
  66. #define MINDOUBLE    DBL_MIN
  67.  
  68. #define MAXF  1.797693134862316
  69. #define MINF  2.225073858507201
  70. #define MAXE  308
  71. #define MINE  (-308)
  72.  
  73. /* another alias for ieee double */
  74. struct ldouble {
  75.     unsigned long hi, lo;
  76. };
  77.  
  78. static int __ten_mul __PROTO((double *acc, int digit));
  79. static double __adjust __PROTO((double *acc, int dexp, int sign));
  80. static double __ten_pow __PROTO((double r, int e));
  81.  
  82. /*
  83.  * mul 64 bit accumulator by 10 and add digit
  84.  * algorithm:
  85.  *    10x = 2( 4x + x ) == ( x<<2 + x) << 1
  86.  *    result = 10x + digit
  87.  */
  88. static int __ten_mul(acc, digit)
  89. double *acc;
  90. int digit;
  91. {
  92.     register unsigned long d0, d1, d2, d3;
  93.     register          short i;
  94.     
  95.     d2 = d0 = ((struct ldouble *)acc)->hi;
  96.     d3 = d1 = ((struct ldouble *)acc)->lo;
  97.  
  98.     /* check possibility of overflow */
  99. /*    if( (d0 & 0x0FFF0000L) >= 0x0ccc0000L ) */
  100. /*    if( (d0 & 0x70000000L) != 0 ) */
  101.     if( (d0 & 0xF0000000L) != 0 )
  102.     /* report overflow somehow */
  103.     return 1;
  104.     
  105.     /* 10acc == 2(4acc + acc) */
  106.     for(i = 0; i < 2; i++)
  107.     {  /* 4acc == ((acc) << 2) */
  108.     asm volatile("    lsll    #1,%1;
  109.              roxll    #1,%0"    /* shift L 64 bit acc 1bit */
  110.         : "=d" (d0), "=d" (d1)
  111.         : "0"  (d0), "1"  (d1) );
  112.     }
  113.  
  114.     /* 4acc + acc */
  115.     asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  116.     asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  117.  
  118.     /* (4acc + acc) << 1 */
  119.     asm volatile("  lsll    #1,%1;
  120.              roxll   #1,%0"    /* shift L 64 bit acc 1bit */
  121.         : "=d" (d0), "=d" (d1)
  122.         : "0"  (d0), "1"  (d1) );
  123.  
  124.     /* add in digit */
  125.     d2 = 0;
  126.     d3 = digit;
  127.     asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  128.     asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  129.  
  130.  
  131.     /* stuff result back into acc */
  132.     ((struct ldouble *)acc)->hi = d0;
  133.     ((struct ldouble *)acc)->lo = d1;
  134.  
  135.     return 0;    /* no overflow */
  136. }
  137.  
  138. #include "flonum.h"
  139.  
  140. static double __adjust(acc, dexp, sign)
  141. double *acc;    /* the 64 bit accumulator */
  142. int     dexp;    /* decimal exponent       */
  143. int    sign;    /* sign flag          */
  144. {
  145.     register unsigned long  d0, d1, d2, d3;
  146.     register          short i;
  147.     register           short bexp = 0; /* binary exponent */
  148.     unsigned short tmp[4];
  149.     double r;
  150.  
  151. #ifdef __STDC__
  152.     double __normdf( double d, int exp, int sign, int rbits );
  153.     double ldexp(double d, int exp);
  154. #else
  155.     extern double __normdf();
  156.     extern double ldexp();
  157. #endif    
  158.     d0 = ((struct ldouble *)acc)->hi;
  159.     d1 = ((struct ldouble *)acc)->lo;
  160.     while(dexp != 0)
  161.     {    /* something to do */
  162.     if(dexp > 0)
  163.     { /* need to scale up by mul */
  164.         while(d0 > 429496729 ) /* 2**31 / 5 */
  165.         {    /* possibility of overflow, div/2 */
  166.         asm volatile(" lsrl    #1,%1;
  167.                     roxrl    #1,%0"
  168.             : "=d" (d1), "=d" (d0)
  169.             : "0"  (d1), "1"  (d0));
  170.         bexp++;
  171.         }
  172.         /* acc * 10 == 2(4acc + acc) */
  173.         d2 = d0;
  174.         d3 = d1;
  175.         for(i = 0; i < 2; i++)
  176.         {  /* 4acc == ((acc) << 2) */
  177.         asm volatile("    lsll    #1,%1;
  178.                  roxll    #1,%0"    /* shift L 64 bit acc 1bit */
  179.                  : "=d" (d0), "=d" (d1)
  180.                  : "0"  (d0), "1"  (d1) );
  181.         }
  182.  
  183.         /* 4acc + acc */
  184.         asm volatile(" addl    %2,%0" : "=d" (d1) : "0" (d1), "d" (d3));
  185.         asm volatile(" addxl   %2,%0" : "=d" (d0) : "0" (d0), "d" (d2));
  186.  
  187.         /* (4acc + acc) << 1 */
  188.         bexp++;    /* increment exponent to effectively acc * 10 */
  189.         dexp--;
  190.     }
  191.     else /* (dexp < 0) */
  192.     { /* scale down by 10 */
  193.         while((d0 & 0xC0000000L) == 0)
  194.         {    /* preserve prec by ensuring upper bits are set before div */
  195.         asm volatile(" lsll  #1,%1;
  196.                     roxll #1,%0" /* shift L to move bits up */
  197.             : "=d" (d0), "=d" (d1)
  198.             : "0"  (d0), "1"  (d1) );
  199.         bexp--;    /* compensate for the shift */
  200.         }
  201.         /* acc/10 = acc/5/2 */
  202.         *((unsigned long *)&tmp[0]) = d0;
  203.         *((unsigned long *)&tmp[2]) = d1;
  204.         d2 = (unsigned long)tmp[0];
  205.         asm volatile(" divu #5,%0" : "=d" (d2) : "0" (d2));
  206.         tmp[0] = (unsigned short)d2;    /* the quotient only */
  207.         for(i = 1; i < 4; i++)
  208.         {
  209.         d2 = (d2 & 0xFFFF0000L) | (unsigned long)tmp[i]; /* REM|next */
  210.         asm volatile(" divu #5,%0" : "=d" (d2) : "0" (d2));
  211.         tmp[i] = (unsigned short)d2;
  212.         }
  213.         d0 = *((unsigned long *)&tmp[0]);
  214.         d1 = *((unsigned long *)&tmp[2]);
  215.         /* acc/2 */
  216.         bexp--;    /* div/2 taken care of by decrementing binary exp */
  217.         dexp++;
  218.     }
  219.     }
  220.     
  221.     /* stuff the result back into acc */
  222.     ((struct ldouble *)acc)->hi = d0;
  223.     ((struct ldouble *)acc)->lo = d1;
  224.  
  225.     /* normalize it */
  226.     r = __normdf( *acc, ((0x03ff - 1) + 53), (sign)? -1 : 0, 0 );
  227.     /* now shove in the binary exponent */
  228.     return ldexp(r, bexp);
  229. }
  230.  
  231. /* flags */
  232. #define SIGN    0x01
  233. #define ESIGN    0x02
  234. #define DECP    0x04
  235. #define CONVF    0x08
  236.  
  237. double strtod (s, endptr)
  238. const register char *s;
  239. char **endptr;
  240. {
  241.     double         accum = 0.0;
  242.     register short flags = 0;
  243.     register short texp  = 0;
  244.     register short e     = 0;
  245.     double       zero = 0.0;
  246.     const char        *tmp;
  247.  
  248.     assert ((s != NULL));
  249.  
  250.     if(endptr != NULL) *endptr = (char *)s;
  251.     while(Isspace(*s)) s++;
  252.     if(*s == '\0')
  253.     {    /* just leading spaces */
  254.     return zero;
  255.     }
  256.  
  257.     if(Issign(*s))
  258.     {
  259.     if(*s == '-') flags = SIGN;
  260.     if(*++s == '\0')
  261.     {   /* "+|-" : should be an error ? */
  262.         return zero;
  263.     }
  264.     }
  265.  
  266.     do {
  267.     if (Isdigit(*s))
  268.     {
  269.         flags |= CONVF;
  270.         if( __ten_mul(&accum, Val(*s)) ) texp++;
  271.         if(flags & DECP) texp--;
  272.     }
  273.     else if(*s == '.')
  274.     {
  275.         if (flags & DECP)  /* second decimal point */
  276.         break;
  277.         flags |= DECP;
  278.     }
  279.     else
  280.         break;
  281.     s++;
  282.     } while (1);
  283.  
  284.     if(Ise(*s))
  285.     {
  286.     if(*++s != '\0') /* skip e|E|d|D */
  287.     {  /* ! ([s]xxx[.[yyy]]e)  */
  288.          tmp = s;
  289.         while(Isspace(*s)) s++; /* Ansi allows spaces after e */
  290.         if(*s != '\0')
  291.         { /*  ! ([s]xxx[.[yyy]]e[space])  */
  292.  
  293.         if(Issign(*s))
  294.             if(*s++ == '-') flags |= ESIGN;
  295.  
  296.         if(*s != '\0')
  297.         { /*  ! ([s]xxx[.[yyy]]e[s])  -- error?? */
  298.  
  299.             for(; Isdigit(*s); s++)
  300.             e = (((e << 2) + e) << 1) + Val(*s);
  301.  
  302.             if(IsValidTrail(*s)) s++;
  303.         
  304.             /* dont care what comes after this */
  305.             if(flags & ESIGN)
  306.             texp -= e;
  307.             else
  308.             texp += e;
  309.         }
  310.         }
  311.         if(s == tmp) s++;    /* back up pointer for a trailing e|E|d|D */
  312.     }
  313.     }
  314.  
  315.     if((endptr != NULL) && (flags & CONVF))
  316.     *endptr = (char *) s;
  317.     if(accum == zero) return zero;
  318.  
  319.     return __adjust(&accum, (int)texp, (int)(flags & SIGN));
  320. }
  321.  
  322. double atof(s)
  323. const char *s;
  324. {
  325. #ifdef __OLD__
  326.     extern double strtod();
  327. #endif
  328.     return strtod(s, (char **)NULL);
  329. }
  330.  
  331.  
  332. /* old naive coding */
  333. #ifdef __OLD__
  334. static double __ten_pow(r, e)
  335. double r;
  336. register int e;
  337. {
  338.     if(e < 0)
  339.     for(; e < 0; e++) r /= 10.0;
  340.     else
  341.     for(; e > 0; --e) r *= 10.0;
  342.     return r;
  343. }
  344.  
  345. #define RET(X)     {goto ret;}
  346.  
  347. double strtod (s, endptr)
  348. const register char *s;
  349. char **endptr;
  350. {
  351.     double f = 0.1, r = 0.0, accum = 0.0;
  352.     register int  e = 0, esign = 0, sign = 0;
  353.     register int texp = 0, countexp;
  354.     
  355.     assert ((s != NULL));
  356.     
  357.     while(Isspace(*s)) s++;
  358.     if(*s == '\0') RET(r);    /* just spaces */
  359.     
  360.     if(Issign(*s))
  361.     {
  362.     if(*s == '-') sign = 1;
  363.     s++;
  364.     if(*s == '\0') RET(r); /* "+|-" : should be an error ? */
  365.     }
  366.     countexp = 0;
  367.     while(Isdigit(*s))
  368.     {
  369.     if(!countexp && (*s != '0')) countexp = 1;
  370.     accum = (accum * 10.0) + Val(*s);
  371.     /* should check for overflow here somehow */
  372.     s++; 
  373.     if(countexp) texp++;
  374.     }
  375.     r = (sign ? (-accum) : accum);
  376.     if(*s == '\0') RET(r);  /* [s]xxx */
  377.     
  378.     countexp = (texp == 0);
  379.     if(*s == '.')
  380.     {
  381.     s++;
  382.     if(*s == '\0') RET(r); /* [s]xxx. */
  383.     if(!Ise(*s))
  384.     {
  385.         while(Isdigit(*s))
  386.         {
  387.         if(countexp && (*s == '0')) --texp;
  388.         if(countexp && (*s != '0')) countexp = 0;
  389.         accum = accum + (Val(*s) * f);
  390.         f = f / 10.0;
  391.         /* overflow (w + f) ? */
  392.         s++;
  393.         }
  394.         r = (sign ? (-accum) : accum);
  395.         if(*s == '\0') RET(r); /* [s]xxx.yyy  */
  396.     }
  397.     }
  398.     if(!Ise(*s)) RET(r);    /* [s]xxx[.[yyy]]  */
  399.     
  400.     s++; /* skip e */
  401.     if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e  */
  402.     
  403.     while(Isspace(*s)) s++;
  404.     if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e[space]  */
  405.     
  406.     if(Issign(*s))
  407.     {
  408.     if(*s == '-') esign = 1;
  409.     s++;
  410.     if(*s == '\0') RET(r); /* [s]xxx[.[yyy]]e[s]  */
  411.     }
  412.     
  413.     while(Isdigit(*s))
  414.     {
  415.     e = (e * 10) + Val(*s);
  416.     s++;
  417.     }
  418.     /* dont care what comes after this */
  419.     if(esign) e = (-e);
  420.     texp += e;
  421.     
  422.     /* overflow ? */ /* FIXME */
  423.     if( texp > MAXE)
  424.     {
  425.     if( ! ((texp == (MAXE+1)) && (accum <= MAXF)))
  426.     {
  427.         errno = ERANGE;
  428.         r = ((sign) ? -HUGE_VAL : HUGE_VAL);
  429.         RET(r);
  430.     }
  431.     }
  432.     
  433.     /* underflow -- in reality occurs before this */ /* FIXME */
  434.     if(texp < MINE)
  435.     {
  436.     errno = ERANGE;
  437.     r = ((sign) ? -HUGE_VAL : HUGE_VAL);
  438.     RET(r);
  439.     }
  440.     r = __ten_pow(r, e);
  441.     /* fall through */
  442.     
  443.     /* all returns end up here, with return value in r */
  444.   ret:
  445.     if(endptr != NULL)
  446.     *endptr = s;
  447.     return r;
  448. }
  449. #endif /* __OLD__ */
  450.  
  451. #else    __M68881__ || sfp004
  452.  
  453. #if 0
  454. #ifndef    sfp004
  455. # define _M68881    /* use the predefined inline functions                */
  456. #endif    sfp004
  457. #endif 0
  458.  
  459. /* M.R's kludgy atof --- 881 version.                        */
  460. /*     uses long integer accumulators and extended precision to put them    */
  461. /*    together in the fpu. The conversion long to extended is done completely    */
  462. /*    on the 881.                                */
  463. /*    using extended precision in _float_ avoids rounding errors.        */
  464.  
  465. /* 12.7.1989, 11.10.90, 28.1.91, 24.11.91                    */
  466. /* On overflow, only +-infinity is returned (the 68881's default).        */
  467. /* 24.11.91: return +-MAXDOUBLE instead of +- INFINITY or NAN            */
  468. /*           set errno to ERANGE/EDOM                        */
  469.  
  470. # include <ctype.h>
  471. # include <stdio.h>
  472. # include <float.h>
  473. # include <math.h>
  474. # include <errno.h>
  475. # include "flonum.h"
  476.  
  477. double atof( const char * );
  478. double strtod( const char *, const char ** );
  479. double _Float_( long, long, long, long );
  480.  
  481. # define    true 1
  482. # define false 0
  483. # define CharIsDigit ( isdigit(*Text) )
  484. # define Digit ((*Text-'0'))
  485.  
  486. #if 0
  487. static unsigned long
  488.     __notanumber[2] = { 0x7fffffffL, 0xffffffffL }; /* ieee NAN */
  489. # define NAN  (*((double *)&__notanumber[0]))
  490. # endif 0
  491.  
  492. # define ten_mul(X)    ((((X) << 2) + (X)) << 1)
  493.  
  494. double strtod( const char * Save, const char ** Endptr )
  495. {
  496.   register int Count; int Negative = false, ExpNegative = false;
  497.  
  498.   double Value;
  499.   union double_di * l_Value;
  500.   register long Exponent, Exp_Temp;
  501.   register long Value_1, Value_2;
  502.   register char c;
  503.   register char * Text;
  504.   register char * Places;
  505.   char Buffer[15];
  506.  
  507.   l_Value = (union double_di *) &Value;
  508.  
  509.   Text = Save;
  510.   Places = Buffer;
  511.  
  512.   /* skip over leading whitespace */
  513.   while (isspace(*Text)) Text++;
  514.  
  515.   if (*Text == '-') {
  516.     Negative = true;
  517.     Text++;
  518.   } else
  519.   if (*Text == '+') {
  520.     Negative = false;
  521.     Text++;
  522.   } else
  523.   if( *Text == 0 ) {
  524.     if( Endptr != NULL ) *Endptr = Text;
  525.     return 0.0;
  526.   }
  527.  
  528.   /* Process the 'f'-part                             */
  529.   /* ignore any digit beyond the 15th                         */
  530.   /* BUG: may overflow if more than 10 digits precede the '.'            */
  531.   /* to cure use to accumulators as is being done for the digits after        */
  532.   /* the '.'.                                    */
  533.  
  534.   Exp_Temp = 0;            /* needed later on for the exponential part    */
  535.   Value_1 = 0; Value_2 = 0; Count = 0; Exponent = 0;
  536.   while( CharIsDigit ) {    /* process digits before '.'            */
  537.     if( Count < 15 ) {
  538.       Count++;
  539.       *Places++ = Digit;
  540.     }
  541.     Text++;
  542.   }
  543.   if ( *Text == '.') {
  544.     Text++;
  545.     while( CharIsDigit ) {    /* process digits after '.'            */
  546.       if( Count < 15 )    {
  547.         Count++;
  548.             *Places++ = Digit;
  549.         Exponent--;
  550.       }
  551.       Text++;
  552.     }
  553.   }
  554.   Places = Buffer;
  555.  
  556.   /* Now, Places points to a vector of <= 15 digits                */
  557.   /* text points to the position immediately after the end of the mantissa    */
  558.   /* Value_2 will contain the equiv. of the 8 least significant digits, while    */
  559.   /* Value_1 will contain the equiv. of the 7 most significant digits (if any)    */
  560.   /* and therefore has to be multiplied by 10^8                    */
  561.   /* no overflow possible in the temporary buffers                */
  562.  
  563.   while( Count > 8 )    {
  564.       Value_1 = ten_mul( Value_1 );      Value_1 += *Places++;
  565.       Count--;
  566.   }
  567.   while( Count > 0 )    {
  568.       Value_2 = ten_mul( Value_2 );      Value_2 += *Places++;
  569.       Count--;
  570.   }
  571.  
  572.   /* 'e'-Part */
  573.   if ( *Text == 'e' || *Text == 'E' || *Text == 'd' || *Text == 'D' ) {
  574.  
  575.     char * Tail = Text;
  576.     Text++;
  577.  
  578.     /* skip over whitespace since ANSI allows space after e|E|d|D        */
  579.     while (isspace(*Text)) Text++;
  580.  
  581.     if ( * Text == '-' ) {
  582.         ExpNegative = true;
  583.         Text++;
  584.     } else
  585.     if( * Text == '+' ) {
  586.         ExpNegative = false;
  587.         Text++;
  588.     }
  589.     if( !CharIsDigit ) {
  590.         *Endptr = Tail;    /* take the 'e|E|d|D' as part of the characters    */
  591.         goto Ende;        /* following the number            */
  592.     } else {
  593.     /* Exponents may have at most 3 digits, everything beyond this will be    */
  594.     /* silently ignored                            */
  595.         Count = 0;
  596.         while( CharIsDigit && (Count < 3) ) {
  597.             Exp_Temp = ten_mul( Exp_Temp ); Exp_Temp += Digit;
  598.               Count++;
  599.             Text++;
  600.         }
  601.         if( ExpNegative ) Exp_Temp = -Exp_Temp;
  602.         Exponent += Exp_Temp;
  603.      }
  604.   }
  605.   Value = _Float_( Value_1, Exponent+8L, Value_2, Exponent );
  606.  
  607.   if( Endptr != NULL ) *Endptr = Text;
  608.   if( (l_Value -> i[0]) >= 0x7ff00000U )    { /* INFINITY or NAN        */
  609.     fprintf(stderr," strtod: OVERFLOW error\n");
  610.     errno = ERANGE;
  611.     Value = DBL_MAX;
  612.   }
  613.  
  614. Ende:    
  615.   if( Negative ) {
  616.     Value = -Value;
  617.   }
  618.   return( Value );
  619.  
  620. # if 0
  621. Error:
  622.   fputs("\njunk number \"",stderr); fputs(Save,stderr); 
  623.   fputs("\" --- returning NAN\n",stderr);
  624.   errno = ERANGE;
  625.   if( Endptr != NULL ) *Endptr = Text;
  626.   return(NAN);    /* == Not A Number (NAN) */
  627. # endif
  628. }
  629.  
  630. double atof( const char * Text )
  631. {
  632.     return(strtod(Text,(char **)NULL));
  633. }
  634.  
  635. /*
  636.  * double _Float_( long Value_1, long Exponent_1, long Value_2, long Exponent_2 )
  637.  *
  638.  * merges the accumulators Value_1, Value_2 and the Exponent to a double
  639.  * precision float
  640.  * called by strtod()
  641.  *
  642.  * does all floating point computations with extended precision on the fpu
  643.  */
  644.  
  645. #endif __M68881__ || sfp004
  646.  
  647. #ifdef __M68881__
  648. asm(
  649. ".even
  650. .text
  651.     .globl __Float_
  652. __Float_:
  653.     ftentoxl    a7@(8),fp1        | load Exponent_1
  654.     ftentoxl    a7@(16),fp0        | load Exponent_2
  655.     fmull        a7@(12),fp0        | fmull Value_2 -> fp0
  656.     fmull        a7@(4),fp1        | fmull Value_1 -> fp1
  657.     faddx        fp0,fp1
  658.     fmoved        fp1,a7@-        | fmoved fp1 -> d0/d1
  659.     moveml        a7@+,d0-d1
  660.      rts
  661. ");
  662. #endif    __M68881
  663. #ifdef    sfp004
  664.  
  665. asm("| mjr, 30.1.1991
  666. |
  667. | base =    0xfffa50
  668. |      the fpu addresses are taken relativ to 'base':
  669. |
  670. | a0: fpu base address
  671. |
  672.  
  673. | waiting loop ...
  674. |
  675. | wait:
  676. | ww:    cmpiw    #0x8900,a0@(resp)
  677. |     beq    ww
  678. | is coded directly by
  679. |    .long    0x0c688900, 0xfff067f8 (a0)
  680. | and
  681. | www:    tst.w    a0@(resp)
  682. |    bmi.b    www
  683. | is coded by
  684. |    .word    0x4a68,0xfff0,0x6bfa        | test
  685. |
  686.  
  687. comm =     -6
  688. resp =    -16
  689. zahl =      0
  690.  
  691.     .globl __Float_
  692. .even
  693. .text
  694. __Float_:
  695.     lea    0xfffa50,a0            | fpu address
  696.  
  697.     movew    #0x4092,a0@(comm)        | ftentoxl -> fp1
  698.     .long    0x0c688900, 0xfff067f8
  699.     movel    a7@(8),a0@            | load Exponent_1
  700.  
  701.     movew    #0x4112,a0@(comm)        | ftentoxl -> fp2
  702.     .long    0x0c688900, 0xfff067f8
  703.     movel    a7@(16),a0@            | load Exponent_2
  704.  
  705.     movew    #0x4123,a0@(comm)        | fmull Value_2 -> fp2
  706.     .long    0x0c688900, 0xfff067f8
  707.     movel    a7@(12),a0@            | load Value_2
  708.  
  709.     movew    #0x40a3,a0@(comm)        | fmull Value_1 -> fp1
  710.     .long    0x0c688900, 0xfff067f8
  711.     movel    a7@(4),a0@            | load Value_1
  712.  
  713.     movew    #0x08a2,a0@(comm)        | faddx fp2 -> fp1
  714.     .word    0x4a68,0xfff0,0x6bfa        | test
  715.  
  716.     movew    #0x7480,a0@(comm)        | fmoved fp1 -> d0/d1
  717.     .long    0x0c688900, 0xfff067f8
  718.     movel    a0@,d0
  719.     movel    a0@,d1
  720.      rts
  721. ");
  722. #endif    sfp004
  723.  
  724. #ifdef TEST
  725. #if 0
  726. #ifdef __MSHORT__
  727. #error "please run this test in 32 bit int mode"
  728. #endif
  729. #endif
  730.  
  731. #define NTEST 10000L
  732.  
  733. #ifdef __MSHORT__
  734. # define    RAND_MAX    (0x7FFF)    /* maximum value from rand() */
  735. #else
  736. # define    RAND_MAX    (0x7FFFFFFFL)    /* maximum value from rand() */
  737. #endif
  738.  
  739. main()
  740. {
  741.     
  742.     double expected, result, e, max_abs_err;
  743.     char buf[128];
  744.     register long i, errs;
  745.     register long s;
  746. #ifdef __STDC__
  747.     double atof(const char *);
  748.     int rand(void);
  749. #else
  750.     extern double atof();
  751.     extern int rand();
  752. #endif
  753.  
  754. #if 0
  755.     expected = atof("3.14159265358979e23");
  756.     expected = atof("3.141");
  757.     expected = atof(".31415"); 
  758.     printf("%f\n\n", expected);
  759.     expected = atof("3.1415"); 
  760.     printf("%f\n\n", expected);
  761.     expected = atof("31.415"); 
  762.     printf("%f\n\n", expected);
  763.     expected = atof("314.15"); 
  764.     printf("%f\n\n", expected);
  765.  
  766.     expected = atof(".31415"); 
  767.     printf("%f\n\n", expected);
  768.     expected = atof(".031415"); 
  769.     printf("%f\n\n", expected);
  770.     expected = atof(".0031415"); 
  771.     printf("%f\n\n", expected);
  772.     expected = atof(".00031415"); 
  773.     printf("%f\n\n", expected);
  774.     expected = atof(".000031415"); 
  775.     printf("%f\n\n", expected);
  776.  
  777.     expected = atof("-3.1415e-9"); 
  778.     printf("%20.15e\n\n", expected);
  779.  
  780.     expected = atof("+3.1415e+009"); 
  781.     printf("%20.15e\n\n", expected);
  782. #endif
  783.  
  784.     expected = atof("+3.123456789123456789"); 
  785.     printf("%30.25e\n\n", expected);
  786.  
  787.     expected = atof(".000003123456789123456789"); 
  788.     printf("%30.25e\n\n", expected);
  789.  
  790.     expected = atof("3.1234567891234567890000000000"); 
  791.     printf("%30.25e\n\n", expected);
  792.  
  793.     expected = atof("9.22337999999999999999999999999999999999999999"); 
  794.     printf("%47.45e\n\n", expected);
  795.  
  796.     expected = atof("1.0000000000000000000"); 
  797.     printf("%25.19e\n\n", expected);
  798.  
  799.     expected = atof("1.00000000000000000000"); 
  800.     printf("%26.20e\n\n", expected);
  801.  
  802.     expected = atof("1.000000000000000000000"); 
  803.     printf("%27.21e\n\n", expected);
  804.  
  805.     expected = atof("1.000000000000000000000000"); 
  806.     printf("%30.24e\n\n", expected);
  807.  
  808.  
  809. #if 0
  810.     expected = atof("1.7e+308");
  811.     if(errno != 0)
  812.     {
  813.     printf("%d\n", errno);
  814.     }
  815.     else    printf("1.7e308 OK %g\n", expected);
  816.     expected = atof("1.797693e308");    /* anything gt looses */
  817.     if(errno != 0)
  818.     {
  819.     printf("%d\n", errno);
  820.     }
  821.     else    printf("Max OK %g\n", expected);
  822.     expected = atof("2.225073858507201E-307");
  823.     if(errno != 0)
  824.     {
  825.     printf("%d\n", errno, expected);
  826.     }
  827.     else    printf("Min OK %g\n", expected);
  828. #endif
  829.     
  830.     max_abs_err = 0.0;
  831.     for(errs = 0, i = 0; i < NTEST; i++)
  832.     {
  833.     expected = (double)(s = rand()) / (double)rand();
  834.     if(s > (RAND_MAX >> 1)) expected = -expected;
  835.     sprintf(buf, "%.14e", expected);
  836.     result = atof(buf);
  837.     e = (expected == 0.0) ? result : (result - expected)/expected;
  838.     if(e < 0) e = (-e);
  839.     if(e > 1.0e-6) 
  840.     {
  841.         errs++; printf("%.14e %s %.14e (%.14e)\n", expected, buf, result, e);
  842.     }
  843.     if (e > max_abs_err) max_abs_err = e;
  844.     }
  845.     printf("%ld Error(s), Max abs err %.14e\n", errs, max_abs_err);
  846. }
  847. #endif /* TEST */
  848.